Probability-generating functions
t0=Sys.time()
#Definition of functions
h <- function(z) {
(1 - z)^((1 - z) / z)
}
px <- function(z,Am,m1,Amm) {
(h(z)^{Am}) * (h(h(z)^m1)^{Amm})
}
psb <- function(z,pb) {
(1 - pb) + pb * z
}
ps <- function(z,p) {
(1 - p) + p * z
}
pgb <- function(z,Am,m1,pb,Amm) {
px(psb(z,pb),Am,m1,Amm)
}
pgs <- function(z,p,Am,m1,Amm){
px(ps(z,p),Am,m1,Amm)
}
#Requiered library for 3D plots
library(plotly)
#Print twenty decimal places
options("digits" = 20)
#Data set to work with
FILENAME="R-GM - FinalPop 2^33, bottle neck 2^20, growth cycles 3, WT2R 10^-7, WT2M 1e-06 M2R 5e-05.txt"
#FILENAME="R-QM - FinalPop 2^33, bottle neck 2^20, growth cycles 3, WT2R 10^-7, WT2M 1e-05 M2R 5e-04.txt"
#FILENAME <- "luria_4trans_1.00e-007mua_2seed.txt"
#Reading the file
td <- read.table(FILENAME,header = T,sep="")
#Sampling probability
ps1 <- 2^{20}/2^{33}
#Dilution process: Binomial sampling, with probability ps1, of resistant cell in the j-th growth cycle
binomial_sampling_GC_j<-function(j,ps1){
a=c()
for (i in 1:nrow(td)){
a=c(a,rbinom(n = 1, size = (td[, (4 * (j - 1) + 2)]+td[ ,(4 * (j - 1) + 4)])[i], prob = ps1))
}
return(a)
}
#Recovery of the k-th probability state of the pmf
#WT2R Mutation rate from wildtype to resistant cells
#WT2M Mutation rate from wildtype to mutator cells
#pop Final population size
#N Total number of probability states
pmf_coeff<-function(k,WT2R,WT2M,pop,N){
sum=0
for(n in 0:(N-1)){
#Auxiliary variable to transform the pgf into a characteristic function
aux_char=exp((2i*pi*n)/N)
#Auxiliary variable to apply the FFT algorithm
aux=(-2i*pi*n*k)/N
#FFT #Fixed scaling factor: 500
sum=sum + pgs(aux_char,p=ps1,Am=10^{WT2R}*pop,m1=10^{WT2M},Amm=10^{WT2R}*pop*500)*exp(aux)
}
#Real part
return(Re(sum/N))
}
#Binomial sampling of resistant cells in the j-th growth cycle, j=1
data<-binomial_sampling_GC_j(j=1,ps1)
#Maximum value in data
N=max(data)
N
## [1] 33
steps=30
#log-likelihood function
ll_fun<-function(WT2R,WT2M){
sum(log(pmf_coeff(data, WT2R,WT2M,pop=2^33,N)))
}
paramters_WT2R=c(-9,-5)
paramters_WT2M=c(-8,-4)
output=c()
for(i in 0:steps){
for(j in 0:steps){
m_WT2R=paramters_WT2R[1]*(paramters_WT2R[2]/paramters_WT2R[1])^{i/steps}
m_WT2M=paramters_WT2M[1]*(paramters_WT2M[2]/paramters_WT2M[1])^{j/steps}
ll=ll_fun(m_WT2R,m_WT2M)
output=rbind(output, c(10^{m_WT2R}, 10^{m_WT2M}, 500*10^{m_WT2R}, ll))
}
}
colnames(output)=c('WT2R', 'WT2M', 'M2R', 'log-likelihood')
head(output)
## WT2R WT2M
## [1,] 1.0000000000000000622e-09 1.0000000000000000209e-08
## [2,] 1.0000000000000000622e-09 1.5230713629662131529e-08
## [3,] 1.0000000000000000622e-09 2.2975616244345304987e-08
## [4,] 1.0000000000000000622e-09 3.4334922203825658472e-08
## [5,] 1.0000000000000000622e-09 5.0841703260677619694e-08
## [6,] 1.0000000000000000622e-09 7.4612270297717860124e-08
## M2R log-likelihood
## [1,] 5.0000000000000008327e-07 -744.49903306033786521
## [2,] 5.0000000000000008327e-07 -744.49115671835045305
## [3,] 5.0000000000000008327e-07 -744.47966875932934272
## [4,] 5.0000000000000008327e-07 -744.46306921362315734
## [5,] 5.0000000000000008327e-07 -744.43930306396680407
## [6,] 5.0000000000000008327e-07 -744.40558319719457359
#Maximum
maximum_GC1=output[order(output[,4], decreasing=TRUE),][1,]
maximum_GC1
## WT2R WT2M
## 7.6856069542669374823e-08 8.0631442982151470284e-05
## M2R log-likelihood
## 3.8428034771334688524e-05 -2.6653791929043785559e+02
#options(width = 10000, height=1000 )
discretization<-function(range, n){
mu_min=range[1]
mu_max=range[2]
a=seq(0,n)
result=mu_min*(mu_max/mu_min)^{a/n}
return(result)
}
paramters_WT2R=c(-9,-5)
WT2R=discretization(paramters_WT2R,steps)
paramters_WT2M=c(-8,-4)
WT2M=discretization(paramters_WT2M,steps)
zz=matrix(output[,4],nrow = steps+1,ncol=steps+1)
x.lab="Mutation rate from WT to R on log<sub>10</sub> scale"
y.lab="Mutation rate from WT to M on log<sub>10</sub> scale"
z.lab="log-likelihood"
main=""
x=-7
y=-6
z=ll_fun(-7,-6)
df_real<-data.frame(x, y, z)
x=log10(maximum_GC1[1])
y=log10(maximum_GC1[2])
z=maximum_GC1[4]
df <- data.frame(x, y, z)
#plot_ly(showscale = T, name=" ", showlegend = TRUE)
fig=plot_ly(x=~WT2R, y=~WT2M, z=~zz, type = 'surface',name= "Log-likelihood surface", showlegend = TRUE, colorbar = list(title = "Log-likelihood"), contours = list(z = list(show=F, usecolormap=TRUE,highlightcolor="#ff0000",project=list(z=T)))) %>% layout(autosize = F, width = 900, height = 650,
title = paste("\n",main),
scene = list(
xaxis = list(title = x.lab),
yaxis = list(title = y.lab),
zaxis = list(title = z.lab)
))%>%add_trace(data = df, x = x, y = y, z = z, mode = "markers", type = "scatter3d",name= "Maximum", showlegend=T ,marker = list(size = 6, color = "red", symbol = 104))
fig=fig%>%add_trace(data = df_real, x = x, y = y, z = z, mode = "markers", type = "scatter3d",name= "Real values", showlegend=T ,marker = list(size = 5, color = "deepskyblue", symbol = 104))
fig
df=as.data.frame(output)
df_confidence<-df[df$`log-likelihood` >= maximum_GC1[4]-qchisq(0.95,2)/2, ]
min(df_confidence$WT2M)
## [1] 1.0000000000000000209e-08
max(df_confidence$WT2M)
## [1] 0.00010000000000000000479
fig1=fig%>%add_surface(showscale = F,name="95% confidence region" ,
contours = list(
z = list(
show = TRUE,
project=list(z=T), # (don't) project contour lines to underlying plane
usecolormap = F, # (don't) use surface color scale for contours
color = "red", # set contour color
width = 30, # set contour thickness
highlightcolor = "black", # highlight contour on hover
start = min(df_confidence$`log-likelihood`), # include contours from z = 0...
end = max(df_confidence$`log-likelihood`), # to z = 1400...
size = max(df_confidence$`log-likelihood`)-min(df$`log-likelihood`) # every 100 units
)
)
)
fig1=fig1%>%layout(autosize = F, width = 900, height = 650)
fig1
#fig %>% add_trace(z = matrix(min(df$`log-likelihood`), nrow = nrow(zz), ncol = ncol(zz)), contours = list(coloring = "none", color = "black"), showscale = F, name="Lower limit of the 95% confidence region")
#fig
#fig%>%add_trace(x=~log10(df$WT2R), y=~log10(df$WT2M), z=~df$`log-likelihood`, type='scatter3d', mode = "markers", marker = list(size = 2))
#fig
# Calculate the posterior probabilities
#df$posterior_prob <- (df$`log-likelihood`) / sum((df$`log-likelihood`))
# Sort the data frame by posterior probability in descending order
#data_frame <- df[order(-df$posterior_prob), ]
# Calculate the cumulative posterior probability
#data_frame$cumulative_prob <- cumsum(data_frame$posterior_prob)
# Identify the values within the 95% credible region
#credible_region <- data_frame[data_frame$cumulative_prob <= 0.975, ]
#credible_region <- credible_region[credible_region$cumulative_prob >= 0.025,]
#fig%>%add_trace(x=~log10(credible_region$WT2R), y=~log10(credible_region$WT2M), z=~credible_region$`log-likelihood`, type='scatter3d', mode = "markers", marker = list(size = 2))
#fig
#options(width = 10000, height=1000 )
discretization<-function(range, n){
mu_min=range[1]
mu_max=range[2]
a=seq(0,n)
result=mu_min*(mu_max/mu_min)^{a/n}
return(result)
}
paramters_WT2R=c(-9,-5)
WT2R=discretization(paramters_WT2R,steps)
paramters_WT2M=c(-8,-4)
WT2M=discretization(paramters_WT2M,steps)
zz=matrix(exp(output[,4]-max(output[,4])),nrow = steps+1,ncol=steps+1)
x.lab="Mutation rate from WT to R on log<sub>10</sub> scale"
y.lab="Mutation rate from WT to M on log<sub>10</sub> scale"
z.lab="Rescaled likelihood"
main=""
x=log10(maximum_GC1[1])
y=log10(maximum_GC1[2])
z=exp(maximum_GC1[4]-max(output[,4]))
df <- data.frame(x, y, z)
#fig=plot_ly(showscale = T, name=" ", showlegend = TRUE)
fig=plot_ly(x=~WT2R, y=~WT2M, z=~zz, type = 'surface',name= "Rescaled likelihood surface", showlegend = TRUE, colorbar = list(title = "Rescaled likelihood"), contours = list(z = list(show=F, usecolormap=TRUE,highlightcolor="#ff0000",project=list(z=T)))) %>% layout(autosize = F, width = 900, height = 650,
title = paste("\n",main),
scene = list(
xaxis = list(title = x.lab),
yaxis = list(title = y.lab),
zaxis = list(title = z.lab)
))%>%add_trace(data = df, x = x, y = y, z = z, mode = "markers", type = "scatter3d",name= "Maximum",
marker = list(size = 6, color = "red", symbol = 104))
fig=fig%>%add_trace(data = df_real, x = x, y = y, z = exp(z-z), mode = "markers", type = "scatter3d",name= "Real values", showlegend=T ,marker = list(size = 5, color = "deepskyblue", symbol = 104))
fig=fig%>%add_surface(showscale = F,name="95% confidence region" ,
contours = list(
z = list(
show = TRUE,
project=list(z=T), # (don't) project contour lines to underlying plane
usecolormap = F, # (don't) use surface color scale for contours
color = "red", # set contour color
width = 20, # set contour thickness
highlightcolor = "black", # highlight contour on hover
start = exp(min(df_confidence$`log-likelihood`)-max(df_confidence$`log-likelihood`)), # include contours from z = 0...
end = exp(max(df_confidence$`log-likelihood`)- max(df_confidence$`log-likelihood`)), # to z = 1400...
size = exp(max(df_confidence$`log-likelihood`)- max(df_confidence$`log-likelihood`)) -exp(min(df_confidence$`log-likelihood`)-max(df_confidence$`log-likelihood`)) # every size units
)
)
)
fig=fig%>%layout(autosize = F, width = 1100, height = 650)
fig